function pli_BL_TASK
inDir='C:\Users\matthew\Desktop\2012\projects\RealAndImaginery\data\';
outDir='C:\Users\matthew\Desktop\2012\projects\RealAndImaginery\results\pli_REST_TASK\';

 batchPLI(inDir,outDir);
function batchPLI(inDir,outDir)


conds={'*_Rleft*','*_Rright*','*_Ileft*','*_Iright*'};
bands={[4 8],[8 13],[13 30],[30 50]};
bandLabels={'theta','alpha','beta','gamma'};
%'*_Rleft*' [4 8], 'theta',
for j=1:length(bands) 
    for i=1:length(conds)
        [restPLI_subjXrepXchXch,taskPLI_subjXrepXchXch]=pliXband(inDir,conds{i},bands{j});
        outfile=strcat(bandLabels{j},conds{i}(2:end-1));
        save(fullfile(outDir,outfile),'restPLI_subjXrepXchXch','taskPLI_subjXrepXchXch');
    end
end


function [restPLI_subjXrepXchXch,taskPLI_subjXrepXchXch]=pliXband(inDir,cond,band)

filter=fullfile(inDir,cond);
files=dir(filter);
nSubjs=length(files);
load(fullfile(inDir,files(1).name),'data');
nCh=length(data.label);
restPLI_subjXrepXchXch=cell(1,nSubjs);
taskPLI_subjXrepXchXch=cell(1,nSubjs);
for i=1:nSubjs
    load(fullfile(inDir,files(i).name),'data');
    bf_data=selectBand(data,band);
    [restPLI_subjXrepXchXch{i},taskPLI_subjXrepXchXch{i}]=pliXtrial(bf_data);
end

function [bf_data]=selectBand(data,band)

bf_data=data;
nTrials=length(bf_data.trial);
for i=1:nTrials
    bf_data.trial{i}=eegfilt(bf_data.trial{i},bf_data.fsample,band(1),band(end));
    bf_data.trialRest{i}=eegfilt(bf_data.trialRest{i},bf_data.fsample,band(1),band(end));
end

function [syncRest,syncTask]=pliXtrial(data)

[nCh,nSample]=size(data.trial{1});
nTrials=length(data.trial);
syncRest=zeros(nTrials,nCh,nCh);
syncTask=zeros(nTrials,nCh,nCh);
%normPLI=zeros(nTrials,nCh,nCh);
for i=1:nTrials
      syncRest(i,:,:)=pli(data.trialRest{i}');
      syncTask(i,:,:)=pli(data.trial{i}');  
end
%normPLI=syncTask./syncRest;
%pliM=squeeze(mean(normPLI,1));

%pliRest=squeeze(mean(syncRest,1));
%pliTask=squeeze(mean(syncTask,1));


function [PLI]=pli(m)

% a is a filtered multichannel signal (time x channels)
% hilbert(a) calculates analytic signal (complex valued) of each
% column of a. Phase Lag Index between channel i and j averaged over
% time bins is stored in PLI(i,j)
% number of channels

nChs=size(m,2);
timeStop=size(m,1);

PLI=zeros(nChs,nChs);
complex_m=hilbert(m);

phi=atan2(imag(complex_m),real(complex_m));


zscoreMatrix=zeros(nChs,nChs);

for i=1:nChs
    
    for j=1:nChs
        count=0;
        count2=0;

        for t=1:timeStop
        
            phaseDiff=phi(t,i)-phi(t,j);
            
            if(sin(phaseDiff)>0)
                count=count+1;
            end
            if(abs(sin(phaseDiff))>0)
                count2=count2+1;
            end
            
        end
        
        if(count2>0)
            pli=2*abs(.5-1*count/count2);
            %for FDR
  %          zscoreMatrix(i,j)=((count/count2)-0.5)/sqrt((0.5*0.5)/count2);
        else
            pli=0;
        end
        
        PLI(i,j)=pli;
    end
end
%sPLI=fdr(zscoreMatrix,PLI,0.01);
function [pliM]=fdr(zscoreMatrix,pliM,alpha)

lindx=find(triu(ones(size(pliM)),1));
[row,col]=ind2sub(size(pliM),lindx);
z=abs(zscoreMatrix(lindx));
k=ones(1,length((z)));%zeros or ones?? stam one
N=length(z);

for i=1:length(z)
    
    k(i)=k(i)+sum(z(i)<z);%to be equal of stam code
    
    if(z2p(z(i))>(k(i)*alpha)/N)
        pliM(row(i),col(i))=0;
        pliM(col(i),row(i))=0;
    end        
        
end
function [p]=z2p(zscore)
        p=1;
       
        if(zscore > 0.67 )
            p=0.5;
        end
        if(zscore > 0.84)
            p=0.4;
        end
        if(zscore > 1.04 )
            p=0.3;
        end
        if(zscore > 1.28 )
            p=0.2;
        end
        if(zscore > 1.64)
            p=0.1;
        end
        if(zscore > 1.96 )
            p=0.05;
        end
        if(zscore > 2.33 )
            p=0.02;
        end
        if(zscore > 2.58 )
            p=0.01;
        end
        if(zscore > 2.81 )
            p= 0.005;
        end
        if(zscore > 3.09 )
            p=0.002;
        end
        if(zscore > 3.29 )
            p=0.001;
        end
        if(zscore > 3.89 )
            p=0.0001 ;
        end


